function [ B_s ] = B_s_fun_4(data_trans_norm,m_regr,combos_all_unique,N,J)

B_s = single(zeros(N,size(combos_all_unique,1)));

for i=1:size(combos_all_unique,1)
     
    ind_other = perms(1:J-1);
        
    combos_x_other = combos_all_unique(i,2:J);
    combos_x_other_r = combos_x_other(ind_other);
    
    combos_z_other = combos_all_unique(i,J+2:2*J);
    combos_z_other_r = combos_z_other(ind_other);
    
    v = single([repmat(combos_all_unique(i,1),[size(ind_other,1),1])   combos_x_other_r...
                repmat(combos_all_unique(i,J+1),[size(ind_other,1),1]) combos_z_other_r]) - 1;
   
    v = unique(v,'rows');       % KEY to avoid counting multiple times the combos that correspond to only one of combos_all        
            
    %v1(i,:) = v; 
     
%     max_combos = 5e6/(N*J^2);
%     thresh = 1:max_combos:size(v,1);
%     thresh = [thresh size(v,1)];
%     
%     temp_B_s = zeros(N,length(thresh)-1);
    
%     for j=1:length(thresh)-1
%         
%         temp_B_s(:,j) = sum(prod(BernsteinElement_2(v(thresh(j):thresh(j+1),:),m_regr,...
%                             [data_trans_norm.x_own data_trans_norm.x_other data_trans_norm.z_own data_trans_norm.z_other]'),3)',2);
%         
%     end
    
    B_s(:,i) = single(sum(prod(BernsteinElement_2(v,m_regr,...
                     [data_trans_norm.x_own data_trans_norm.x_other data_trans_norm.z_own data_trans_norm.z_other]'),3)',2));

end
 
end